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ABSTRACT 


A  method  is  devised  for  the  transient-state  thermal 
stress  analysis  of  two  pipes,  joined  butt  to  butt  and 
subjected  to  rapid  or  sudden  change  of  internal  fluid 
temperature.   Although  it  is  assumed  that  there  is 
symmetry  about  the  common  pipe  axis,  the  properties  of 
the  materials  as  well  as  the  thickness  of  each  pipe  may 
be  different.   Then,  given  a  specified  time-temperature 
relationship  for  the  internal  fluid,  over  a  specified 
problem  time,  the  maximum  stress  encountered  may  be 
obtained.   A  digital  computer  program  is  appended  for 
the  solution  of  such  problems. 
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I  .   INTRODUCTION 

In  the  design  and  analysis  of  certain  piping  of  critical 
importance,  such  as  nuclear  power  piping,  approved  methods, 
such  as  those  in  the  USAS  B31.7  Code,  call  for  analysis  of 
transient  thermal  stresses.   The  most  important  cause  of 
such  transient  thermal  stresses  is  variation  of  the 
temperature  of  the  fluid  flowing  in  the  piping.   It  is 
conventional  to  regard  the  fluid  as  having  infinite  heat 
capacity  and  to  presume  that  the  fluid  flows  through  the 
piping  with  sufficient  rapidity  that  all  parts  of  the  piping 
component  or  assembly  in  question  are  simultaneously 
exposed  to  the  same  interior  ambient  temperature.   At  any 
rate,  an  analysis  based  upon  such  simplifying  assumptions 
may  provide  the  engineering  answer  that  is  required  or  may 
form  the  basis  for  a  more  precise  evaluation. 

The  geometry  of  piping  components  and  assemblies  is 
generally  rather  complicated  and  one  method  of  analysis  is 
to  study  a  greatly  simplified  geometry  and  then  to  apply 
"stress  indices"  which  are  tabulated  constants  intended  to 
reflect  the  greater  severity  of  the  thermal  stresses  in  the 
actual  object  as  compared  to  that  having  the  idealized 
geometry.   One  such  simplified  geometry  of  importance  is 
that  of  two  semi-infinite  pipes,  having  a  common  axis,  and 
joined  together  by  welding  or  brazing  at  the  origin  of  an 
axial  coordinate  system.   If  the  two  pipes  have  identical 


dimensions  and  properties,  determination  of  the  transient 
thermal  stresses  is  a  rather  simple  problem:  there  is  no 
variation  along  the  axial  length.   However,  if  the  two 
pipes  do  not  match  in  thickness  or  in  material  properties, 
so-called  dislocation  stresses  will  be  set  up  near  the 
junction  and  the  task  of  determining  the  significant 
transient  thermal  stresses  becomes  considerably  more 
difficult.   It  may  be  noted  that  the  analysis  to  be  given 
herein  also  contemplates  a  mismatch  in  pipe  diameter,  but 
this  is  not  important  since  the  theory  employed  will  be 
limited  to  the  case  of  thin-wall  pipes. 

The  problems  of  interest  arise  upon  considering  the 
ambient  fluid  temperature,  to  which  the  interior  of  the 
piping  is  convectively  exposed,  as  varying  with  respect 
to  time  in  a  specified  fashion.   In  the  development  which 
follows,  the  general  approach  is  capable  of  dealing  with 
any  temperature-time  relation.   The  implementation  presented 
is  in  the  form  of  an  immediately  useful  digital  computer 
program  for  the  analysis  of  these  problems.   In  addition, 
there  is  specialization  to  temperature-time  relations 
which  postulate  a  long  period  at  constant  temperature 
followed  by  a  sudden  change  in  temperature  or  by  a  rapid 
linear  variation  in  temperature  which  is  in  turn  followed 
by  maintenance  of  a  new  steady  temperature. 

The  analysis  of  the  thermal  stresses  in  any  case  con- 
sists of  the  following  steps: 
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(1)  determination  of  the  temperature  distribution 
in  the  piping 

(2)  determination  of  the  thermal  stresses  in  the  piping 
at  a  considerable  distance  from  the  junction 

(3)  determination  of  the  dislocation  reactions,  i.e., 
the  radial  force  and  bending  moment  (per  unit 
circumference)  exerted  by  one  of  the  semi-infinite 
pipes  upon  the  other,  and  the  calculation  of 
stresses  resulting  therefrom 

(4)  combining  these  stresses  in  an  appropriate  manner 
to  find  what  could  be  called  the  "significant 
stress"  and  locating  the  point  in  the  piping  and 
the  point  in  time  for  which  this  significant  stress 
is  greatest 

The  task  is  greatly  simplified  because  the  problem  at 
hand  is  one  possessing  symmetry  about  the  common  pipe  axis. 
It  is  assumed  that  the  fluid  temperature  is  the  same  func- 
tion of  time  at  all  points  adjacent  to  the  inner  wall  of 
the  pipe  and  that  the  coefficient  of  convective  heat 
transfer  is  likewise  the  same  at  all  such  points.   It  is 
also  assumed  that  the  exterior  surfaces  of  both  pipes  are 
perfectly  insulated  and  that  the  following  properties  of 
both  pipes  are  given  and  classified  into  two  categories: 
first,  thermal  conductivity,  specific  heat,  and  density 
which  may  vary  with  temperature;  and  second,  Young's 
modulus  of  elasticity,  Poisson's  ratio,  and  the  coefficient 
of  linear  thermal  expansion  which  must  remain  constant. 
Under  certain  conditons  the  properties  in  the  second  cate- 
gory might  be  permitted  to  vary,  but  this  has  not  been 
investigated  in  this  thesis.   It  is  also  necessary  to  know 
pipe  wall  thickness  and  outside  diameter. 
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The  task  of  determining  the  transient  temperature 
distribution  in  each  of  the  two  pipes  is  simplified  by 
the  assumption  that  there  is  no  axial  flow  of  heat  even 
though  adjacent  points  in  the  two  pipes  may  not  be  at  the 
same  temperature.   In  other  words,  an  error  is  introduced 
by  neglecting  axial  flow  of  heat  near  the  junction  but  the 
problem  is  too  difficult  to  handle  without  this  assumption. 
There  is  no  evident  way  to  estimate  the  effect  of  the 
assumption  but  the  effect  should  be  small  for  thin  wall 
pipes  to  which  this  entire  development  is  aimed. 

With  this  assumption,  the  determination  is  relatively 
straightforward.   If  the  properties  of  the  materials  were 
not  themselves  dependent  upon  temperature,  there  might  be 
some  incentive  to  seek  analytic  solutions.   However,  for 
the  sake  of  generality,  a  numerical  method  has  been  used. 
This  consists  of  slight  modifications  to  a  digital  computer 
program,  BETTY1,  developed  by  John  E.  Brock,  Professor  of 
Mechanical  Engineering,  Naval  Postgraduate  School,  Monterey, 
California.   This  program  is  based  upon  the  general  algo- 
rithm described  in  Reference  5,  but  including  several 
auxiliary  features;  see  Appendix  B. 

The  output  of  BETTY1  is  a  tabulation,  at  each  instant 
in  time,  of  the  temperature  at  equally  spaced  radii  through 
the  pipe  wall.   In  the  present  application  this  provides 
too  much  data  to  handle  so  that  the  next  step  is  to  pass  a 
polynomial  approximation  of  chosen  degree  (usually  of  third 
or  fourth  degree)  through  the  data  generated  by  BETTY1 ,  and 
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retain  the  coefficients  of  this  polynomial  for  each  instant 
of  time.  This  also  provides  for  more  convenient  processing 
at  later  stages  of  the  analysis. 

Next,  using  this  information,  the  stress  distribution 
and  radial  displacement  is  calculated  for  each  pipe  at  a 
position  remote  from  the  junction.   This  is  a  simple  and 
standard  sort  of  calculation  which  is  simplified  by  virtue 
of  the  temperature  data  being  given  in  polynomial  form. 
Here,  and  elsewhere  in  what  follows  such  calculations  must 
be  made  for  each  pipe  at  each  instant  in  time. 

The  next  step  is  the  determination  of  the  dislocation 
reactions,  i.e.,  the  radial  force  and  bending  moment  per 
unit  circumference  at  the  pipe  junction.   This  calculation 
is  based  upon  the  fact  that  continuity  of  radial  deforma- 
tion and  its  first  derivative  with  respect  to  the  axial 
coordinate  must  be  maintained  at  the  junction.   It  leads  to 
a  pair  of  simultaneous  equations  in  which  the  results  of 
the  calculation  described  in  the  preceding  paragraph  enter. 
Having  the  dislocation  reactions,  the  corresponding  disloca- 
tion stresses  may  be  determined.   These  all  contain  a  nega- 
tive exponential  term  and  die  out  away  from  the  junction. 
However,  they  are  not  necessarily  maximum  right  at  the 
junction,  so  that  some  care  is  required  later  to  be  sure 
that  the  worst  condition  is  not  overlooked.   The  actual 
stress  components  at  any  point  are  obtained  by  superposing 
the  thermal  stresses  remote  from  the  junction  and  the  dis- 
location stresses  near  the  junction. 
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Finally  it  is  necessary  to  determine  the  "significant 
stress"  or  measure  of  stress  severity.   In  conformity  with 
the  stress  criteria  used  by  the  most  applicable  codes,  the 
maximum  shearing  stress  criterion  has  been  employed  here. 
At  any  cross  section,  at  any  radius,  and  at  any  time,  this 
maximum  shearing  stress   may  be  calculated  without  great 
difficulty.   Because  the  variation  with  respect  to  radius 
occurs  only  in  the  non-dislocation  stresses  (except  for  a 
linear  variation)  and  because  the  variation  with  respect  to 
axial  position  occurs  only  in  the  dislocation  stresses,  the 
problem  of  determining  a  maximum  shearing  stress  is 
simplified. 

In  brief,  the  analysis  and  its  implementation  in  the 
form  of  a  digital  computer  program  called  PAM1  perform  all 
the  above  steps  and  emerge  with  a  statement  of  the  magnitude 
of  the  "most  significant  stress,"  its  location,  and  the 
time  at  which  it  occurs. 

The  effort  in  developing  this  thesis  has  been  in  under- 
standing and  relating  the  several  steps  as  described  above, 
in  writing  sub-programs  to  perform  the  operations  in  each 
of  the  steps  as  outlined  above,  making  use  of  programs 
(such  as  BETTY1  and  a  least  squares  curve  fitting  routine) 
that  appeared  to  be  of  immediate  use,  debugging  and  proving 
each  portion,  knitting  the  separate  parts  together  into  a 


It  is  conventional  to  calculate  and  employ  a  quantity 
known  as  the  "significant  stress"  or  the  "stress  intensity" 
which  is  equal  to  twice  the  maximum  shearing  stress.   This 
practice  is  adhered  to  in  what  follows  herein. 
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complete  whole,  and  finally  debugging  and  proving  the  final 
program. 

The  steps  of  the  theoretical  analysis  and  the  methods 
devised  to  employ  them  are  described  in  Sections  II  and  III 
hereof.   The  digital  computer  program  is  described,  in 
general  terms,  in  Section  IV.   Instructions  for  the  use 
of  the  program  are  given  in  Appendix  A.   The  program  list- 
ing is  given  in  Appendix  B. 

There  was  not  sufficient  time  after  developing  the 
program  to  make  a  systematic  use  of  it  to  investigate  the 
effect  of  certain  material  or  geometrical  parameters  upon 
the  value  of  the  significant  stress.   It  is  recommended 
that  this  become  the  point  of  departure  for  a  future  thesis 
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II.   THEORY 

When  the  wall  of  a  pipe  is  nonuniformly  heated  as  for 
instance  when  it  is  subjected  to  a  rapid  or  sudden  change 
of  internal  fluid  temperature,  its  elements  do  not  expand 
uniformly  and  mutual  interference  sets  up  thermal  stresses. 
It  is  assumed  that  the  distribution  of  the  temperature  does 
not  vary  along  the  axis  of  the  pipe  and  that  there  is 
complete  symmetry  about  the  axis.   A  ring  is  cut  from  the 
pipe  by  two  plane  cross  sections  perpendicular  to  the  axis 
and  a  unit  distance  apart.   During  thermal  deformation  such 
cross  sections  can  be  assumed  to  remain  plane  if  taken 
sufficiently  distant  from  the  ends  of  the  pipe.   Therefore, 
the  unit  elongation  in  the  direction  of  the  axis  is  constant 
[7].   Then,  the  unit  elongations  in  the  perpendicular 
directions  are 


dv; 
e   =  -=—  =  constant 
z    dz 


z      =—  (1) 

er    dr 


t    r 


These  elongations  can  be  expressed  as  functions  of  the 
axial,  radial,  and  circumferential  stresses  remote  from  the 
end  of  the  pipe  (o  ,  a  ,  and  o      respectively) ,  the  thermal 
expansion,  a,  and  the  temperature,  T,  where  T  varies  with 
the  radial  distance,  r,  only. 
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From  Reference  6,  p.  66,  the  elongations  are 


a 

e„  =  -# £-  (ar  +  a,  )  +  aT 

z  F     r     t 


£r  "  ^  "  ~t    (oz  +  °t>  +  aT  <2> 


ct  -  %  -  -g-  (0Z  +  or)  +  aT 


Let  A  be  the  increase  in  unit  volume.   From  this 


A  =  ez  +  er  +  et  =  1  "E2p  (az  +  ar  +  afc)  +  3  aT      (3) 


Substituting  equations  (2)  into  equation  (3)  it  can  be 
shov/n  that 


E   /         y   a \     cTE 


z    1+y    z     1  -  2y  '    1  -  2y 


»r  ■  i§u  (£r +  r^>  -  i^ir 


le.  +   ^-A>  -   aTE 


t    1  +  y    t     1  -  2y      1  -  2y 

From  Peference  7,  p.  206,  the  equation  for  radial  equi' 
librium  of  an  element  is 

da     a     -  a. 

r     r    t 

Substituting  equations  (4)  and  (1)  into  equation  (5) 

2 
dulduul+ydT  ,,.. 

dr2  r      dr  2         1  -  y         dr 
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which  may  be  written  in  the  form 


dF  [r   d7  (ru)]  =  I" ^  a  ^  (7) 


Integration  with  respect  to  r  gives 


~  (ru)  =  }  + y  aTr  +  2C,r  (8) 

dr         1  -  y  1 


A  second  integration  gives 


2T 

u  =  -i-  y4-J|  /   aTr  dr  +  C,r  +  C~  -i-  (9) 

a. 


where  C,  and  C?  are  constants  of  integration.   C,  and  C~    are 
found  from  the  boundary  conditions  at  the  inner  and  outer 
surfaces  of  the  pipe,  namely 


(a  )     =  0  ;   (a  )     =  0  (10) 

r  r=a  r  r=b 


From  equations  (10)  the  constants  of  integration  are  eval- 
uated as 

2     b 

C   =  t      -=-^ 5-  /   aTr  dr 

2    1  -  JJ  ,  2     2 
z  b   -  a   a 


_  (1  +  y)  (1  -  2y)     1     r      nTr   Hr  _  „_         n-n 

Cl (1  -  y) 72 2  /  a   r   r    y£z        (11) 

D   —  a   a 


An  expression  for  a   is  obtained  by  substituting  e  =  du/dr 

and  e   =  u/r  into  the  second  equations  (4)  and  then  taking 

u  from  equation  (9)  and  the  values  of  the  constants  from 
equations  (11) .   This  gives 
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r   ,   r  2     2b      1 

/   aTr  dr  +  -£ =-^ *-  /  aTr  dr  <       (12) 


r    1  -  y       2  2 „  2     2X 

r   a  r  (b   -  a  )  a        I 


An  expression  for  a   is  obtained  from  the  equation  of 
equilibrium  (5),  thus 


d 
o^   =   o      +   r  -t—  a 
t     r     dr   r 


From  this , 

r  2     2b  i 


E 
a,  = 


t    1  -  y 


—z     f      aTr  dr  +   ?r   *  a  -   /  aTr  dr-  aT 
r2   a  r  (b^  -  a^)  a 


(13) 


An  expression  for  a   is  obtained  by  substituting  £   =  du/dr 
and  e.  =  u/r  into  the  first  of  equations  (4)  and  then  taking 
u  from  equation  (9)  and  the  values  of  the  constants  from 
equation  (11).   This  expression  will  contain  the  constant 
elongation,  e  ,  in  the  direction  of  the  axis  of  the  pipe. 
If  it  is  assumed  that  the  pipe  can  expand  freely  in  the 
axial  direction,  the  magnitude  of  £   can  be  calculated 
from  the  condition  that  the  sum  of  the  axial  forces  over 
any  normal  cross  section  is  equal  to  zero.   From  this 

2       b 
ez  =  -z /   aTr  dr  (14) 

b   -  az      a 


and 

»2  =  - fh 27   '    aTr  dr  -  r^n  (15) 

(1  -  y)  (b   -  a  )   a 
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Substituting  the  values  of  the  constants  from  equation  (11) , 
the  radial  displacement  remote  from  the  pipe  junction  is 

2 


u  =-yrez  +  (r^rr) 


,   r        (1  -  2y)r  +—  b 

-  f    aTrdr  + ^ 2 — ~  f    aTrdr  !   (16) 

a  b   -  a     a       : 


Thus  the  thermal  stresses,  cj  .  o    .    and  o        and  radial 

r    t       z 

displacement,  u ,  at  a  distance  from  the  pipe  junction  may 
be  determined  by  use  of  equations  (12),  (13),  (15)  and  (16). 
It  should  be  noted  here  that  at  a  distance  remote  from  the 
pipe  junction  there  is  no  variation  of  the  stresses  or 
radial  displacement  with  respect  to  z,  the  axial  coordinate 
of  the  pipe.   There  exists  here  a  bending  moment  per  unit 
circumference,  denoted  as  M^,  where 

b  ,-V 

MM=/   a_  (r  -  2_  +  _D)  dr  (17) 

a  c. 

See  Figure  1 ,  case  A 

On  the  other  hand,  if  the  end  of  the  pipe  is  supported 
or  clamped,  as  for  example  when  it  is  joined  to  another 
pipe,  free  expansion  of  the  pipe  is  prevented  and  a  bending 
moment  per  unit  circumference,  M   and  a  radial  force  per 
unit  circumference,  V  ,  are  set  up  at  the  pipe  junction 
(Figure  1,  case  B) .   Because  of  the  continuity  with  the 
joining  pipe,  the  moment,  stresses,  and  radial  deflection 
at  the  pipe  junction  will  be  different  from  those  at  a 
distance  from  the  end  of  the  pipe.   In  other  words,  the 
adjoining  pipe  exerts  a  bending  moment  and  radial  force. 
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It  should  be  noted  here  that  the  total  bending  moment  at 

the  pipe  junction  is  M   -  M  ,  which  will  be  represented 

here  as  M  .   Thus,  in  solving  for  the  dislocation  reactions, 
o 

although  M   is  obtained  directly,  it  is  worth  noting  that 
that  bending  moment  depends  on  both  the  moment  at  the  pipe 
junction  as  well  as  that  remote  from  the  junction.   The 
total  radial  force  is  just  that  at  the  pipe  junction,  V  , 
since  there  is  no  force  remote  from  the  junction  (Figure  1). 
This  bending  moment  and  radial  force  cause  a  deflection 
and  tip  slope  at  the  end  of  the  pipe  (z  =  0) ,  which  are 
denoted  by  n  and  9,  respectively.   The  same  is  true,  however 
for  the  other.   Therefore,  let  n, /  and  6, ,  denote  the  deflec- 
tion and  tip  slope,  respectively,  at  the  end  of  the  adjoin- 
ing pipe.   Now  if  these  pipes  are  to  join  up  continuously, 
the  total  deflection  for  each  pipe  at  the  pipe  junction 
must  be  the  same.   The  total  deflection  consists  of  the 

deflection  due  to  M  and  V   as  well  as  the  radial  displace- 

o      o  r 

ment,  u,  remote  from  the  junction.   Equating  these  total 
deflections,  thus 

n  +  u  =  n-,  +  u±  (18) 

Also,  if  these  pipes  are  to  join  continuously,  the  tip 
slopes  at  the  ends  of  each  pipe  must  be  negative  of  each 
other.   That  is, 


6  =  -  61  (18a) 
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Figure  1.   Combination  of  Cases  A  and  B 
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From  Reference  7,  it  is  seen  that 


-Bz 
n  =  e  5   t$M   (-sin  Bz  +  cos  Bz)  +  V   cos  Bz]       (19) 
26  D    ° 


where 


12(1  -  y2) 


and 


4    3(1-  y2) 

p     ,a+b,2u2 
(— )  h 


and  where  the  algebraic  signs  have  been  changed  to  conform 
to  the  convention  herein  that  n  is  positive  outward.   If 
the  pipe  junction  is  taken  as  the  origin  of  the  axial 
coordinate  (i.e.,  z  =  0) ,  then  the  deflection  at  the  junc- 
tion is 


n  =  -V-   [BM'  +  V  ]  (20) 

2B  D  ° 


Also  from  Reference  1 ,    the  tip  slope  is  the  first  derivative 
of  the  deflection  with  respect  to  z,  or 


dn  ~2z 

~  =  -=-s—    [2^M'    cos    Bz    +   V       (cos    Bz    +    sinBz)]     (21) 

dz     2g2D     o  o 


Therefore  the  tip  slope  at  the  pipe  junction  (z  =  0)  is 

9  = —■  [2BM*  +  V]  (22) 

2B2D     °     ° 
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It  should  be  noted  that  the  radial  displacement,  u, 
given  in  equation  (16)  is  a  function  of  the  radius,  r. 
The  value  of  u  used  in  equation  (18)  was  the  average  value 
of  the  radial  displacement  found  at  the  inner  and  outer 
radii  of  the  pipe.   This  seems  to  be  the  most  reasonable 
and  convenient  of  several  different  choices  that  could  be 
made . 

Noting  that  6  and  D  are  functions  of  the  material  and 

the  geometry  of  each  pipe,  their  values  may  be  obtained  and 

substituted  into  equations  (20)  and  (22).   Also  the  values 

of  u  and  u, ,  may  be  obtained  from  equation  (16).   For 

mechanical  equilibrium,  we  must  have  M  =  M  -  M  =M  ,  =  M  ,-M  , 

o    o   °°   ol  ol   °°1 

and  V  =  -V  , .   These  relations,  together  with  equations 
o      ol  v'" 

(18)  and  (18a)  permit  solving  for  the  dislocation  reactions 

M  ,  M*  ,  and  V  .   Hence,  the  dislocation  stresses  near  the 
o    ol       o 

pipe  junction  can  be  obtained. 

From  Reference  7,  the  axial  stress  due  to  the  bending 
moment  is 

where  y  is  the  radial  distance  measured  from  midsurface,  so 
that  -  -  <  y  <  -.   Also  from  Reference  7 

M(z)=-£  [e"ezsin  Bz]  +  m'  [e"3z  (cos  Bz+sin3z)]   (24) 
p  o 

The  maximum  axial  stress  is  obtained  at  the  axial 
position  where  M(z)  is  a  maximum.   Since  M(z)  contains  a 
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negative  exponential  term,  it  will  eventually  decrease  to 
zero  as  z  increases.   But  this  does  not  necessarily  mean 
that  M(z)  will  be  maximum  at  the  origin.   The  method  of 
finding  the  point  of  maximum  bending  moment  is  shown  in 
detail  in  the  next  section. 

From  Reference  7  the  circumferential  stress  due  to  the 
dislocation  reactions  is 

=  En_(zj_  + 

tz     r        zz 

m 

Where  n(z)  comes  from  equation  (19)  and  r   =  .   It  should 

^  m     2 

be  noted  here  that  there  is  no  radial  stress  due  to  the 
dislocation  reactions  because  of  the  assumption  that  the 
pipe  wall  is  thin. 

Summarizing  now,  there  is  a  radial,  circumferential, 
and  axial  stress  remote  from  the  junction  as  well  as 
additional  circumferential  and  axial  stresses  near  the 
junction.   Each  of  these  is  some  function  of  the  radius, 
so  that  for  a  given  radius  and  given  axial  location  these 
stresses  may  be  superposed  on  each  other  and  the  total 
stresses  found.   Using  equations  (12),  (13),  (15),  (23), 
and  (25)  these  total  stresses  are  defined  as 


r 
°Tfc  =  °t  +  °tZ  (26) 

°T   =  °z  +  °zz 
z 

The  leading  subscript  T  indicates  "total." 
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The  radial  force  V,  also  gives  rise  to  a  shearing 

stress,  t    .   Presuming  that  this  shearing  stress  is 

Trz  3 

distributed  parabolically ,  it  can  be  represented  as 


Trz 


3V  rn    ,2y_  2 


2h 


[1 


h 


(27) 


where,  from  Reference  7 


V  =  e  ez  [V   (cos  3z  -  sin  Bz)  -  2Bm'  sin  3z]      (28) 


Thus  a  general  element  is  loaded  as  shown  in  Figure  2. 
(It  should  be  noted  that  this  element  is  also  subjected  to 
the  circumferential  stress,  a   ) . 


Tr 


Figure  2 
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The  three  principal  stresses  are  then 


al  "  aTt 


QT2+   aTr  /,gTz  "gTrx2  2 

a2  =  = +  /( 5 }   +  TTrz  (29) 


gT2+aTr    ;,qTz  "°Tr,2        2 
a3  =  2 '  ( 2 }   +  TTrZ 


\ 


The  maximum  significant  stress  at  a  point  is  defined  as 


a  .   =  Max  {la,  -  a~ | , I a0  -  cJ  ,|a,  -  a,  I  }  .        (29a) 
sig         '1     2 '  '  2     3 '  '  3     l1 


It  is  necessary  to  find  the  point  in  space  (i.e.,  values  of 

r  and  z)  and  the  point  in  time  for  which  a  .   achieves  its 

^  sig 

maximum  value,  and  to  determine  this  value.   This  must  be 
done  for  each  of  the  two  pipes  individually  since  allowable 
values  of  a  .   may  be  different  for  the  two  pipes.   How  this 
is  done  is  considered  in  detail  in  Section  III. 
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III.   METHOD 

As  indicated  in  Section  I,  the  program  BETTY1  has  been 
incorporated  into  the  program  PAM1  to  obtain  the  transient 
temperature  distribution.   This  is  a  program  for  dealing 
with  the  one  dimensional  thermal  transient  analysis  of 
homogeneous  slabs  or  tubes.   Its  use  herein  is  for  the 
thermal  transient  analysis  of  a  homogeneous  pipe  in  which 
the  outer  surface  is  perfectly  insulated  and  the  inner 
surface  is  convectively  exposed  to  a  fluid  having  a 
specified  time-temperature  relationship.   Physical  prop- 
erties of  the  pipe  material  (i.e.  thermal  conductivity, 
weight  density,  and  specific  heat)  may  be  temperature 
dependent;  the  program  is  set  up  to  consider  property- 
temperature  relationships  approximated  by  a  third  degree 
polynomial.   A  similar  variation  is  permitted  in  the 
surface  heat  transfer  coefficient. 

The  following  is  a  brief  summary  of  the  general 
algorithm  of  the  program  [5].   As  shown  in  Figure  3  a  unit 


Figure  3 
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length  of  pipe,  contained  in  a  small  dihedral  angle,  is 
considered  to  be  divided  into  subvolumes ,  where  the  inner 
and  outer  subvolume  have  half  the  thickness  of  the  others. 
The  thermal  properties  of  each  subvolume  are  considered  to 
be  concentrated  at  its  central  nodal  point,  and  heat  is 
imagined  to  be  conducted  between  nodal  points  through  a 
network  of  fictitious  heat-conducting  rods  of  appropriate 
thermal  conductance.   In  the  transient  state,  however, 
heat  is  not  only  conducted  to  and  from  each  nodal  point  in 
the  lumped  network,  but  additionally  each  nodal  point 
experiences  a  change  in  internal  energy  in  an  amount  which 
depends  upon  its  temperature  change  during  a  given  time 
increment,  its  specific  heat,  the  total  volume  of  material 
which  it  represents,  and  the  density  of  the  material. 

Assuming  that  the  specific  heat,  c,  and  the  thermal 
conductivity,  K,  are  uniform,  the  heat  conducted  by  each 
rod  is  -KAT,  where  AT  is  the  temperature  difference 
between  adjacent  nodal  points.   The  total  heat  conducted 
in  the  finite  time  increment  At  is  -KATAt.   Now  the  change 
in  internal  energy  of  a  given  nodal  point  is  cpVAT  ,  where 
p  is  the  material  density,  V  is  the  subregion  volume,  and 
T   is  the  temperature  of  the  nodal  point  at  the  end  of 
the  finite  time  interval  At.   The  heat  balance  can  be 
written  as 

-ZKATAt  =  cpVAT'  (30) 

Using  an  example  of  three  nodal  points,  as  shown  here 
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X 


C,p,V,K 


Figure  4 


a  specific  example  of  equation  (30)  is 


K(T2  -  T]_)  +  K(T3  -  T1)  =  ^  (TJ  -  T±) 


Solving  for  T' ,  the  future  temperature  of  the  nodal  point,. 


T|  =  y[T2  +  T3  +  (  — 2)  Tx] 


(31) 


where 


Y  = 


KAt 
cpV 


and  is  known  as  the  Fourier  modulus.   From  this  it  can  be 
seen  that  the  future  temperature  of   1   depends  on  its 
present  temperature  and  the  present  temperature  of  its 
adjacent  nodal  points.   Equation  (31)  is  the  basis  of  the 
transient-state  numerical  method.   It  should  be  noted 
that  the  primary  purpose  of  BETTY1  was  to  find  the  quan- 
tities delta-tee-one  and  delta-tee-two,  defined  and  required 
by  the  Nuclear  Power  Piping  code,  but  in  the  present 
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application  only  the  transient-state  temperature  distri- 
bution is  of  interest. 

The  immediately  applicable  accessible  output  from 
BETTY1  is  a  set  of  nodal  point  temperatures  at  each  epoch 
of  time.   For  convenience  in  what  follows,  it  was  decided 
to  represent  this  temperature  distribution  by  a  polynomial 
in  r,  using  least  squares  procedures.   This  results  in  an 
economy  of  computer  storage  and  convenience  in  performing 
the  integrations  required  in  equations  (9)  and  (11) 
(et.seq.),  with  minor  loss  of  accuracy.   A  brief  discussion 
of  the  curve  fitting  procedure  follows. 

The  deviation  at  a  data  point  is  defined  as  the  differ- 
ence between  the  tabulated  T-value  and  the  T-value  computed 
for  a  functional  relationship.  The  least-squares  criterion 
asks  for  a  minimum  of 


N  2 

G  =  Z   [T1  (r.  )  -  T.  ]  (32) 

,11 
1  =  1 


where  N  equals  the  number  of  data  points,  T   denotes  a 
value  computed  at  r  =  r.  from  a  functional  relationship, 
and  T.  denotes  the  corresponding  tabulated  value.   This 
can  be  differentiated  to  find  its  minimum,  which  leads  to 
equations  that  in  many  cases  of  practical  interest  are 
linear  and,  in  principle  at  least,  easy  to  solve. 

The  form  of  the  relationship  must  be  chosen  in  advance 
Since  it  is  important  to  find  a  relationship  which  will  be 
easy  to  obtain  and  integrate,  a  polynomial  of  the  form 
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T(r)  =  a   +  a,r  +  a2r   + 


k    „      i 
.  .  a,  r  =  Z   a  .rJ 
k     .     -j 


(33) 


was  selected.   Here  k  is  the  degree  of  the  polynomial  and 
the  coefficients  a-.,  a,,  a~  ...  a,  must  be  determined  to 
minimize  G. 

Substituting  equation  (33)  into  equation  (32) 


N     k      .       „ 
G  =  Z   (  Z   a.r?  -  T. ) 
i=l   j=0   ^  1 


(34) 


For  a  minimum,  the  partial  derivative  of  G  with  respect 

to  each  coefficient,  a  ,  (m  =  0 ,  1 ,  2  k)  must  vanish. 

m 

Thus 


0  =  |£_  =   z  [2(1      a.r?  -  T. )  rm] 

3a     .  ,      .  A   i  l     l    iJ 
m   i=l      j  =  0   J 


N 


.J+I1^  _ 


N 


m. 


=  2  [  Z   a.  (  Z  r->  '"')  -   Z   T.r  ] 

■  n         3   •   i    !         -i    11 

-)  =  0   J  1=1  1  =  1 


(35) 


m  =  0,  1,  2 
This  set  of  equations  may  be  arranged  as 


N 


Zr.   Zr 
i     l 


2     3 
Zr.  Zr    Zr 
li     l 


Zr . 

i 

v  k+1 

Zr . 

l 


-  k  v  k+1  „  k+2    2k 
Zr  .  Zr .    Zr r . 

ii     l      i 


"                        " 

a 
o 

ZT. 

l 

al 

= 

ZT.  r. 
l    l 

ak 

ZT.r. 
l    l 

i 

.             m 

(36) 
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N 
where,  in  this  last  expression,  Z  means   Z  .   It  is 

i=l 
straightforward  to  calculate  the  coefficients  and  to  solve 

the  system  to  obtain  the  desired  coefficients  a  . 

2  m 

For  the  sudden  or  rapid  temperature  changes  which 
appear  to  be  of  greatest  interest,  it  seems  likely  that 
the  degree  of  the  polynomial  need  not  exceed  k  =  4 ,  al- 
though no  verification  of  this  conjecture  has  been  under- 
taken.  However,  it  should  be  noted  here  that  even  for 
some  unusual  time-temperature  relationship  such  as  a 
temperature  which  fluctuates  with  time,  a  polynomial  of 
larger  degree  could  still  be  used  to  represent  this 
temperature  distribution. 

Now,  for  each  pipe  and  for  each  time  epoch,  the  co- 
efficients of  the  polynomial  representation  of  the 
temperature  as  a  function  of  the  radius  are  available. 
These  polynomials  can  be  integrated  and  the  stresses 
remote  from  the  pipe  junction  can  be  determined.   Equa- 
tion (12) 


ar  =  iqi  I"  4  /   aTrdr  +   2*2"  a  2,  '  aTrdr]      (12) 
M     ra  r  (b  -  a  )  a 


will  serve  as  an  example  here.   Given  values  of  E,  y ,  a,  a, 
and  b,  and  the  particular  radius,  r,  as  well  as  the  poly- 
nomial for  the  temperature  distribution,  the  radial  stress 
can  be  found.   On  a  computer,  the  radius  can  be  incremented 
from  the  inner  radius,  a,  to  the  outer  radius,  b,  thereby 
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finding  the  radial  stress  at  increments  throughout  the 
thickness  of  the  pipe.   A  similar  procedure  is  followed 
using  equations  (13)  and  (15),  thereby  finding  the  circum- 
ferential and  axial  stresses  throughout  the  thickness  of 
the  pipe . 

It  was  noted  earlier  that  the  radial  displacement/ 
remote  from  the  pipe  junction,  used  in  determining  the 
total  deflection  would  be  an  average  value  of  the  radial 
displacements  found  at  the  inner  and  outer  surfaces  of 
the  pipe ,  i.e., 

u  =  |  [u(a)  +  u(b)]  (37 

From  this  it  can  easily  be  shown  that 


u  =  (£±£)  [-ye   +  22(1+  V\    f   aTrdr]  (38) 

b   -  a   a 


where  e   is  found  from  equation  (14).   Thus,  this  radial 
displacement  can  be  determined  for  each  pipe  by  evaluating 
the  integral  of  the  polynomial  for  the  temperature  dis- 
tribution at  the  inner  and  outer  surfaces  of  the  pipe. 

The  next  step  is  to  find  the  bending  moment,  M   and 
the  radial  force  V   at  the  junction  of  the  two  pipes.   This 
procedure  was  outlined  briefly  in  the  preceding  section, 
but  will  be  shown  in  some  detail  here.   Substituting  equa- 
tion (20)  and  the  values  obtained  for  u  and  u,  ,  from 
equation  (38) ,  into  equation  (18)  thus 
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[3M'  +  V  ]  +  u  =  1 [PtM'   +  V  ]  +  u. 


2B3D    °    °        23^    1  ol    o      1 


where  as  before,  subscript  1  is  used  to  distinguish  between 

the  two  pipes.   Now  substituting  M   =  M   +  M   and 
c   c  ^   o    o    °° 


M'   =  M   +  M  , .  and  rearranging 
ol     o     °°1  r   r 


i                i                           11                                      MM, 
[_J_ i_]M     +  t_l i_j  v        u      u_  [.J! |i_j        (39) 

23   D      26p1  23   D      2pp1  23   D      23^ 


Substituting  equations  (22)  into  equation  (18a) 


[23M1  +  V  ]  =  — ^ [23-.M',  +  V  ] 


232D      °     °     20^     1  ol     ° 


then  substituting  as  above  and  rearranging 


M     M 


Equations  (39)  and  (40)  are  easily  solved  for  the  mechanical 
reactions  M   and  V   after  which  the  appropriate  dislocation 
reaction  moments 


M'  =  M   +  M 
o     o     °° 


M  .  =  M   +  M  . 
ol     o     °°1 


may  be  determined. 

Now  that  the  dislocation  reactions,  M  ,  M  , ,  and  V  , 

o    ol        o 

have  been  determined,  the  stresses  which  they  produce  in 
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each  pipe  can  be  computed.   These  stresses  may  then  be 
superposed  on  the  thermal  stresses  already  obtained  and 
shown  in  equations  (12),  (13)  and  (15),  obtaining  the  total 
stresses  as  defined  in  equation  (26).   Next  computing  the 
shearing  stress,  equation  (27),  and  the  principal  stresses, 
equation  (29),  the  significant  stress  may  be  determined. 
It  is  worth  repeating  here  that  this  significant  stress 
is  the  largest  difference  between  any  pair  of  principal 
stresses.   Hence,  it  is  a  function  of  the  radial  coordinate, 
as  well  as  the  axial  coordinate  of  the  pipe.   Also  it  may 
be  useful  to  know  which  pair  of  principal  stresses  produce 
this  maximum. 

The  computer  program  takes  the  absolute  value  of  the 
difference  of  each  pair  of  the  three  principal  stresses, 
and  then  selects  the  largest  of  these  as  the  maximum 
significant  stress.   This  is  done  at  increments  of  the 
radius  throughout  the  thickness  of  each  pipe.   It  has  been 
noted  earlier  that  the  radial  stress  is  zero  at  the  inner 
and  outer  surfaces  of  the  pipe.   Also  the  circumferential 
and  axial  stresses  remote  from  the  pipe  junction  are  maximum 
at  either  the  inner  or  outer  surface  for  a  large  class  of 
fluid  temperature  variations.   Therefore,  it  seems  reason- 
able that  in  taking  differences,  the  largest  difference 
might  be  at  either  inner  or  outer  surface.   In  fact,  in 
certain  important  examples  of  specified  time-temperature 
relationships  (i.e.  sudden  change  in  internal  temperature), 
it  has  been  observed  empirically  by  the  writer  that  this 
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maximum  does  occur  at  either  the  inner  or  outer  surface, 
but  admittedly  this  may  not  always  be  the  case.   In  any 
case,  finding  the  radial  position  at  which  the  significant 
stress  is  maximum  is  simplified  by  merely  finding  the  sig- 
nificant stress  at  given  increments  of  radius  throughout 
the  pipe  thickness. 

The  axial  position  at  which  the  significant  stress  is 
maximum  is  that  position  at  which  the  dislocation  stresses 
are  maximum.   As  noted  earlier,  these  stresses  are  maximum 
at  the  point  of  maximum  bending  moment  which  may  occur  at 
the  junction  (z  =  0)  or  at  a  mathematical  extremum  of  the 
bending  moment.   Equation  (24)  gives  the  moment  as  a  func- 
tion of  z,  3/  M1  and  V  .   But  now  3,  m'  and  V   are  known. 

o      o  O      o 

Therefore,  taking  the  derivative  of  the  moment  with  respect 
to  z  and  setting  it  equal  to  zero,  the  points  of  extremal 
moment  can  be  determined. 

V  _ft 

M(z)    =   -£    [e    pzsinBz]    +   M*     [e~pz (cos    Bz    +    sin3z)]       (24) 

p  o 

Setting  the  derivative  equal  to  zero, 

0  =  Vq  [e"ez(cos  3z  -  sin  3z)]-  3M'[e~ez(2  sin  3z) ] 

Rearranging  to  a  form  readily  usable  on  the  computer,  and 
solving  for  z 

1  Vo 

z  =     Arctan  [    ,°+    ]  (41) 

o    o 
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This  method  is  employed  in  the  computer  program.   In 
this  manner,  the  first  (and  only  physically  meaningful) 
point  of  extremal  moment  can  be  found  for  each  pipe ,  being 
sure  to  use  the  appropriate  value  of  3.   It  must  be 
remembered  here,  that  the  value  of  the  dislocation  moment 
at  the  pipe  junction,  M1  ,  has  already  been  determined.   In 
examining  certain  important  temperature-time  histories, 
such  as  those  corresponding  to  sudden  or  rapid  changes  in 
internal  fluid  temperature,  it  has  been  observed  empirically 
that  the  maximum  bending  moment  for  each  pipe  occurs  at 
the  pipe  junction.   From  these  examples,  it  was  noted 
that  the  magnitude  of  the  moment,  found  by  substituting 
the  value  of  z  obtained  from  equation  (41)  into  equation 
(24),  was  always  less  than  the  magnitude  of  the  moment  at 
the  pipe  junction.   In  fact,  in  all  the  examples  studied 
in  this  thesis  investigation,  not  only  is  this  moment 
at  the  junction  a  maximum,  but  the  moment  decreases 
rapidly  away  from  the  junction,  oscillating  above  and  below 
a  value  of  zero.   The  program  uses  either  the  end  moment 
or  the  moment  at  the  axial  position  given  by  equation  (41) , 
whichever  is  larger,  in  computing  the  dislocation  stresses. 

If  the  temperature  distribution  were  such  that  the 
maximum  significant  stress  did  not  occur  at  the  inner  or 
outer  surfaces,  the  user  of  program  PAM1  still  has  the 
option  of  wholesale  determinations  at  different  radii  from 
which  a  maximum  may  be  found.   If  the  maximum  significant 
stress  did  not  occur  at  the  end  (z  =  0),  the  preceding 
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determination,  which  is  based  upon  maximizing  M  without 
considering  n  or  V  would  probably  not  give  the  axial 
position  of  greatest  significant  stress.   Because  the  cases 
which  seemed  to  be  of  greatest  interest  lead  to  maximum 
significant  stress  occuring  at  either  the  inner  or  outer 
surface,  the  program  PAM1  in  its  present  form  does  not 
include  the  capability  of  making  determinations  at  the 
axial  positions  other  than  for  (z  =  0)  and  the  value  of 
z  given  in  equation  (41).   Accordingly,  for  cases  in  which 
the  maximum  significant  stress  occurs  between  the  inner 
and  outer  radii,  the  program  does  not  presently  provide  a 
capability  of  finding  the  critical  spot  and  determining 
the  true  maximum  significant  stress.   However,  there  should 
be  no  great  difficulty  in  modifying  the  program  to  provide 
this  capability  if  desired. 

Once  the  maximum  significant  stress  is  found  for  a 
given  time,  the  same  procedure  can  be  followed  for  the  next 
time  increment.   It  may  be  worth  repeating  that  two  values 
of  maximum  significant  stress  are  determined  at  each  given 
time  epoch-one  for  each  pipe.   The  temperature  distribution 
is  known  for  that  time  increment,  therefore,  a  polynomial 
can  be  found  to  represent  the  temperature,  and  then  that 
polynomial  can  be  used  to  find  the  next  significant  stress. 
This  method  is  employed  in  the  computer  program  making  it 
possible  to  find  a  maximum  significant  stress  over  any  time 
period. 
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IV.   DIGITAL  COMPUTER  PROGRAM 

This  section  will  serve  as  a  summary  for  all  that  has 
preceded  it  in  the  sense  that  now  all  the  theoretical 
techniques  and  the  methods  devised  to  carry  out  these 
techniques  have  been  incorporated  into  the  program  PAM1 
which  is  immediately  usable  for  the  solution  of  the  partic- 
ular problem  at  hand.   This  program  is  written  in  Fortran  IV 
language  and  employs  only  normal  Fortran  functions. 

Appendices  A  and  B  will  give  a  complete  listing  of  the 
program,  as  well  as  instructions  for  its  use.   The  writer's 
purpose  in  this  section  is  to  outline  the  program  (a  flow 
chart  in  words  perhaps)  thereby  describing  the  various 
sub-programs  employed  and  woven  together  to  form  the  final 
program 

As  indicated  in  the  preceding  section,  the  whole  idea 
of  applying  the  theoretical  analysis  devised  herein  to  a 
transient-state  problem  rests  on  the  premise  that  a  tran- 
sient temperature  distribution  could  be  obtained,  and  it  was 
shown  that  by  using  a  numerical  method,  such  a  distribution 
can  be  obtained.   In  view  of  the  fact  that  this  temperature 
had  to  be  integrated  in  order  to  find  the  stresses  which 
it  produces,  the  temperatures  and  the  corresponding  radii 
found  from  BETTY1 ,  were  used  as  data  points;  then  a  poly- 
nomial was  fitted  to  these  points.   The  program  employs  the 
least  squares  method  and  Gauss  elimination  to  find  the 
coefficients  of  the  polynomials. 

40 


PAM1  reads  data  for  one  pipe  and  solves  for  the  co- 
efficients of  the  polynomial  for  every  instant  in  time  at 

which  a  temperature  distribution  has  been  determined. 

2 

This  continues  until  the  problem-end-time   has  been  reached, 

with  the  program  storing  each  of  the  coefficients  thus  far 
obtained.   Next,  the  program  reads  data  for  the  joining 
pipe  and  solves  for  the  coefficients  of  a  similar  polynomial. 
Now,  however,  at  each  time  interval  the  program  recalls 
the  coefficients  for  the  other  pipe  at  the  corresponding 
time.   The  program  continues  to  do  this  until  the  compu- 
tation terminates.   For  each  time  increment  there  is  a 
polynomial  for  the  temperature  distribution  for  each  pipe; 
hence,  the  stresses  can  be  determined  for  every  epoch  of 
time . 

As  shown  in  Section  II,  these  stresses  are  determined 
by  integrating  the  polynomial,  which  is  a  function  of 
radius,  multiplied  by  the  radius  (i.e.  /Trdr) .   This  process 
is  carried  out  in  subroutine  PINT.   Subroutine  PVAL  is  then 
used  to  evaluate  the  integral  at  any  given  radius,  as  well 
as  at  the  inner  and  outer  surfaces  of  the  pipe.   In  cal- 
culating the  circumferential  and  axial  stresses,  the 
temperature  at  a  given  radius  is  required.   Subrouting  TVAL 
evaluates  only  the  polynomial  for  the  temperature  distribu- 
tion in  order  to  obtain  the  temperature  at  a  given  radius. 


2 

This  is  an  upper  limit,  introduced  as  an  input  to 

PAM1 .   The  computation  terminates  when  time  has  been  incre- 
mented to  the  value  of  problem-end-time. 
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The  same  techniques  are  employed  in  finding  the  radial 
displacement  remote  from  the  pipe  junction.   Subroutines 
PINTl,  PVALl,  and  TVALl  perform  the  same  procedures  in 
finding  the  stresses  and  radial  displacement  for  the 
adjoining  pipe. 

The  program  then  solves  equations  (39)  and  (40)  for 
the  unknowns  M   and  V  ,  the  total  bending  moment  and  radial 
force  at  the  pipe  junction,  and  the  corresponding  disloca- 
tion reactions  (M1 ,  V  )  and  (M',  V  )  are  next  determined. 

o    o        ol    o 

At  this  point,  equation  (41)  is  used  to  determine  the  point 
of  maximum  bending  moment.   The  moment  at  this  point  is 
then  computed  and  compared  to  the  moment  at  the  pipe 
junction,  and  the  one  of  greater  magnitude  is  used  in 
determining  the  dislocation  stresses.   The  total  stresses 
are  then  determined  and  the  maximum  significant  stress 
as  defined  in  Section  II  is  obtained. 

Each  of  the  stresses  obtained  thus  far  depends  upon 
the  radius,  so  that  the  program  can  find  these  stresses 
for  each  pipe  at  any  radius  throughout  the  thickness  of 
the  pipe.   Since  it  has  been  this  writer's  experience 
that  the  maximum  significant  stress  occurs  at  either  the 
inner  or  outer  surface,  the  program  offers  the  option  of 
either  finding  only  those  stresses  at  the  inner  or  outer 
radii  or  all  stresses  at  given  increments  of  radius.   This 
depends  on  the  wish  of  the  user  and  is  elected  through 
the  use  of  a  data  card. 
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In  any  case,  once  these  stresses  are  determined  for  a 
given  time  the  program  continues  to  the  next  instant  in 
time,  calculates  the  coefficients  of  the  polynomial  for  the 
second  pipe,  recalls  the  coefficients  already  determined 
for  the  first  pipe,  and  performs  all  the  evaluations 
described  above. 

These  procedures  continue  until  the  problem-end-time 
is  reached.   After  determining  the  significant  stress  for 
each  epoch  of  time,  the  program  calculates  the  largest  of 
these,  thereby  determining  the  maximum  significant  stress 
encountered  over  the  problem  time  for  each  pipe.   PAM1  also 
furnishes  the  time  and  location  of  this  maximum  stress, 
as  well  as  indicating  the  nature  of  the  two  principal 
stresses  responsible  for  it.   It  should  be  repeated  that 
the  output  includes  two  such  combinations  of  results  -  one 
for  each  pipe. 
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APPENDIX  A 
Instructions  for  use  of  Program  PAM1 

This  program  is  designed  for  the  IBM  series  360  com- 
puter; however,  there  is  very  little  machine  dependence. 
Testing  and  execution  thus  far  has  been  on  an  360/67  con- 
figuration.  All  quantities  beginning  with  the  alphabetic 
letters  I,  J,  L,  M,  and  N  are  integer  variables,  used  for 
distinguishing  specific  locations  in  the  pipe  wall,  used 
as  counters  for  storing  purposes,  and  used  as  limits  for 
DO-loops.   All  other  quantities  are  continuous  (REAL  *8, 
i.e.,  double  precision)  variables.   Many  of  these  variables 
differ  in  appearance  only  by  the  number  (1)  placed  at  the 
end  of  the  variable.   This  is  to  distinguish  related 
variables  for  each  of  the  joining  pipes.   For  instance, 
ELAS  is  used  to  denote  the  modulus  of  elasticity  for  one 
pipe,  while  ELAS1  denotes  the  modulus  of  elasticity  for 
the  joining  pipe. 

The  program  determines  temperatures  at  the  centers  of 
each  of  the  N  layers  into  which  the  pipe  wall  is  divided. 
The  arrangement  of  these  layers  is  one  half-layer,  (N-l) 
full  layers,  and  a  final  half-layer.   An  example  for  N=5 
layers  is  illustrated  in  Figure  A-l. 
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Figure  A-l 

Note  that  the  centers  are  numbered  from  C,  to  C.T ,  ,  .   There- 

1      N+l 

fore,  in  the  program  the  given  value  of  N  is  increased  by 
unity  so  as  to  permit  numbering  the  centers  as  illustrated 
in  the  figure  above.   One  important  consequence  of  this  is 
that  in  fitting  the  polynomial  to  the  temperature  distri- 
bution, there  are  N+l  temperatures  and  N+l  radii  used  as 
data  points.   The  number  of  layers  should  be  between  four 
and  twenty,  inclusive,  but  if  greater  than  ten,  it  should 
be  an  even  number.   Otherwise  the  program  automatically 
changes  N  to  conform  to  these  conditions. 

The  temperatures,  coefficients  of  the  polynomials,  and 
related  stresses  are  calculated  at  discrete  epochs,  uni- 
formly separated  by  time-interval  DT.   The  given  value  of 
DT,  as  selected  by  the  user  is  reduced  if  necessary  to 
assure  convergence  of  the  computation.   As  this  calculated 
value  may  be  different  for  each  of  the  connected  pipes, 
care  must  be  taken  to  see  that  the  smaller  value  is  used. 
Because  of  the  sequence  of  operations  involved,  it  may 
happen  that  temperature  evaluations  have  been  completed 
for  one  pipe  before  analysis  of  the  parameters  relating 
to  the  second  pipe  has  revealed  that  a  shorter  time 
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interval  should  have  been  used.   In  this  case  an  appropriate 
message  is  printed  and  the  computation  ceases;  it  is  then 
necessary  for  the  user  to  make  input  modifications  as 
advised  and  to  restart  the  program.   However,  if  this 
case  does  not  arise,  the  computation  procedes  to  completion 
as  desired.   Another  method  of  insuring  that  the  same  value 
of  DT  is  used  for  both  pipes  is  to  remove  temporarily  the 
card  DT  =  BB  which  precedes  statement  31.   This  enforces 
the  input  value  of  DT.   However,  if  this  is  done,  one 
should  be  alert  to  the  possibility  of  instability  or  non- 
convergence  in  the  calculation. 

The  data  deck  is  arranged   in   two  parts.   The  first 
part  contains  material  and  geometrical  properties  which 
are  constant  for  each  pipe.   This  part  consists  of  two 
data  cards.   The  second  part  of  the  data  deck  contains  the 
numbers  of  layers  into  which  the  pipe  is  divided,  time 
interval,  problem-end-time,  data  for  those  properties  of 
the  material  and  fluid  which  may  vary  with  temperature, 
as  well  as  information  regarding  the  internal  fluid  time- 
temperature  relationships  and  initial  pipe  temperature. 
The  degree  of  the  polynomial  desired  should  also  be 
selected  at  this  time.   Although  there  has  not  been 
sufficient  experimental  use  of  the  program  to  warrant 
making  definite  statements,  tentatively  it  may  be  recom- 
mended that  the  degree  of  the  polynomial  should  not 
exceed  0.6  N  (where  N  is  the  number  of  layers).   Unless 
the  fluid  temperature  fluctuates  with  time  there  seems  to 
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be  no  reason  to  employ  a  polynomial  of  higher  than  fourth 
degree.   Another  important  input  is  the  option  parameter 
NSTAR  which  determines  if  the  stresses  will  be  calculated 
at  only  the  inner  and  outer  surfaces  or  at  given  increments 
of  radius  throughout  the  wall  of  the  pipe.   Presently 
evaluations  are  made  for  r  =  a  +  n(b-a)/100;  n=0 ,  1,...,100. 
However  the  interval  may  be  modified  by  changing  cards  900 
and  901. 

Detailed  instructions  for  the  preparation  of  each  data 
card  follow. 

Part  I  Card  1 

Column  Format  Description  or  Value 

1-10  F  10.0  Young's  modulus  of  elasticity, 

psi 

11-20  F  10.0  coefficient  of  linear, 

thermal  expansion,  °F 

21-30  F  10.0  Poisson's  ratio 

31-40  F  10.0  inner  radius  of  pipe,  inches 

41-50  F  10.0  outer  radius  of  pipe,  inches 

Card  2 

This  card  is  punched  exactly  as  Card  1  using  parameters 
now  for  the  connecting  pipe. 

In  the  following,  note  that  all  integer  inputs  must  be 
right-justified  in  their  respective  fields. 
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Part  II  Card  3 


Column  Format  Description  or  Value 

1-3  13  degree  of  polynomial  desired 

4-6  13  1  or  2.   Number  one  (1)  is 

used  to  denote  first  pipe; 
number  two  (2)  is  used  to 
denote  second  pipe. 

7-9  13  1  or  2 .   Option  parameter  one 

(1)  used  to  calculate 
stresses  at  inner  and  outer 
surface  only,  while  option 
parameter  two  (2)  calculates 
stresses  at  increments  of 
radius  throughout  pipe  wall. 


Card  4 

1  II  1.  (A  digit  in  this  position, 

cards  5  through  11  merely 
sets  up  appropriate  internal 
branching) 

2-3  12  number  of  layers  into  which 

pipe  thickness  is  divided 

4-5  12  1.   Number  one  (1)  must  be 

punched  here.   Program  BETTYl 
had  capabilities  to  print 
output  only  on  certain  time 
intervals.   Program  PAM1  dees 
not. 

6-10  15  blank 

11-20  F10.4  outside  diameter  of  pipe, 

inches 

21-30  F10.4  pipe  thickness,  inches 

31-40  F10.4  time  interval,  seconds 

41-50  F10.4  maximum  problem  time,  seconds 
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Column 

1 

11-20 

Card  5 

Format 

11 

2 

F10.4 

s 

Description  or  Value 


surface  heat  transfer  co- 
efficient, Btu/hr  -  sq.ft  - 
°F 

21-30  F10.4  Higher  coefficients  which 

31-40  give  variation  with  respect 

41-50  to  temperature  according  to 

the  formula;    Property  = 
F5  +  X*  (F6  +  X*  (F7  +X*F8)  ) 
where  X  denotes  temperature 
and  F5  through  F8  denote  the 
four  entries  in  columns  11 
through  50  respectively.   If 
the  higher  coefficients  are 
left  blank,  a  value  of  zero 
is  used. 


Card  6 

1  II  3 

11-20  F10.4  thermal  conductivity,  Btu/hr 

-  ft  -  °F 

21-30  F10.4  higher  coefficients  used  as 

31-40  F10.4  in  Card  5 

41-50  F10.4 


Card  7 

1  II  4 

11-20  F10.4  specific  heat,  Btu/lb  -  °F 

21-30  F10.4  higher  coefficients  used  as 

31-40  F10.4  in  Card  5 

41-50  F10.4 
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Card  8 

Column  Format              Description  or  Value 

1  II  5 

11-20  F10.4  density,  lbs/cu.ft. 

21-30  F10.4  higher  coefficients  used  as 

31-40  F10.4  in  Card  5 

41-50  F10.4 


Card  9 

1  II  6 

11-20  F10.4  time,  seconds,  at  which  the 

internal  fluid  temperature 
is  that  entered  in  columns 
21-30 

21-30  F10.4  temperature  of  internal 

fluid,  °F.   Type  6  cards 
govern  the  time-temperature 
relationship  of  the  internal 
fluid  and  as  many  as  thirty 
may  be  used.  These  can  be 
denoted  as  Cards  9b,  9c, 
etc. 


Card  10 

1  II  7 

11-20  F10.4  Initial  uniform  temperature 

of  pipe,  °F 


Card  11 


II  8   Number  eight  (8)  signals 

computer  to  begin  computation 


Part  II  of  the  data  deck  must  now  be  repeated  using  the 
same  formats  but  this  time  using  those  parameters  for  the 
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connecting  pipe  (i.e.,  the  pipe  identified  in  Card  2). 
These  cards  follow  directly  after  the  card  with  the  number 
eight  (8)  punched  in  column  1,  and  they  must  also  end  with 
a  card  with  8  punched  in  column  1. 
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APPENDIX  B 


Glossary  of  Symbols  Used  in  Program 


Listing  of  program  starts  on  page  59 


A(I) 

A3 

AA,BB3,RR 

AAA 

AENDX 

AENDMO 

AL(I) 

ARG 

ALPHA 

AYE 

B 

B3 

BBB3 

BBX 

BEST 

BETA 

BIG 
BUZZ  (I) 


cross  sectional  area  at  center  number  I 

inner  radius  of  pipe 

quantities  useful  in  calculating  /Trdr 

quantity  useful  in  calculating  SIGMAR 

absolute  value  of  ENDX 

absolute  value  of  ENDMO 

K(X(I)  )/(SPHT(X(I)  )  *DENS(X(I)  )  )  ,  a 
useful  combination 

value  for  which  integral  is  evaluated 

coefficient  of  linear  thermal  expansion 

input  from  data  card 

input  from  data  card 

outer  radius  of  pipe 

quantity  useful  in  calculating  SIGMAR 

value  of  BETA  times  z 

maximum  significant  stress  for  a  given 
time 


3(1  -  y  ) 


2,2 

r  h 


largest  value  of  RT ,  TZ  and  ZR 

variable  used  to  store  radius  of  maximum 
significant  stress  for  each  time  interval 
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BX 

C 

C3 

ccc 

D 

D3 

DENOM 

DENS  (Z) 

DENA, 

DENB, 

DENC, 

DEND 

DEL 

DELTA1 

DELTA2 
n,PT,rn 

JJ£j  1  1  J. 

DETT2 
DONE, 

DTWO 

DD 

DDD 

DR 

DT 

E 

ELAS 

EMINF 

ENDX 

ENDP 


BBX 

input  from  data  card 

value  of  /Trdr  from  inner  radius  to  some 
radius ,  R 

quantity  useful  in  calculating  SIGMAT 

input  from  data  card 

value  of  /Trdr  from  inner  radius  to 
outer  radius 

a  quantity  useful  in  calculating  TAVG 

metal  density  at  temperature  z 

coefficients  for  determining  DENS(Z) 

3        2 
plate  f lexural  rigidity  =  Eh  /12  (1  -  y  ) 

AT,    for   USAS   B31.7   Nuclear  Power  Piping  Code 

AT      for   USAS   B31.7   Nuclear  Power  Piping  Go de 

absolute  value  of  DELTA1 

absolute  value  of  DELTA2 

quantities  used  in  interpolating  in 
subroutine  FLUTEM 

quantity  useful  in  determining  radial 
displacement 

quantity  useful  in  obtaining  SIGMAZ 

thickness  of  a  layer 

time  interval 

input  from  data  card 

modulus  of  elasticity 

bending  moment  remote  from  pipe  junction 

bending  moment  as  a  function  of  z 

dislocation  bending  moment 
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ENDMO 

EPSZ 

F 

F3 

FELTA1 

FELTA2 

FLUI (I) ,  FLUO(I) 

FLUTEM  (Z,A,B) 

FRAC 

FLU 
G 

GAKA 
GAUSS 2 

GLUG1 

GLUG2 

H 

H3 

HO 

IDIMT 
IDIMX 


total  bending  moment  at  pipe  junction 
unit  elongation  in  axial  direction 
input  from  data  card 

quantity  useful  in  calculating  SIGMAR 
cumulatively  largest  value  of  DETT1 

cumulatively  largest  value  of  DETT2 

3 
interior,  exterior  fluid  temperature 

subroutine  entered  with  time  Z  and 
giving  A,B  as  corresponding  inside, 
outside  fluid  temperatures 

numerical  constant  used  in  making 
initial  estimate  of  inner  surface 
temperature 

interior  fluid  temperature 

quantity  useful  in  calculating  SIGMAT 

radius  of  maximum  significant  stress 

subroutine  used  for  finding  coefficients 
of  polynomial  for  temperature  distri- 
bution 

value  of  DETT1  at  time  TT1 

value  of  DETT2  at  time  TT2 

interior  surface  film  coefficient 

quantity  useful  in  calculating  SIGMAZ 

external  surface  film  coefficient  (taken 
as  zero  corresponding  to  perfect 
insulation) 

number  of  coefficients  in  temperature 
polynomial 

number  of  coefficients  in  an  associated 
temperature  polynomial 


Reference  to  exterior  fluid  temperature  is  for  a  planned 
future  capability  of  BETTY1;  no  such  quantity  is  actually 
used  in  the  calculations  herein  described. 
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IDIMY 
JAP 
JUAY 
K(Z) 

KA  ,  KB  ,  KC  ,  KD 

KAY ( I ) 

LORDER 

LL 

MO 

MOM 

MOP 

Nx/N2/N3/N4 
NO 

NCOUNT 

NFLAG 

NPRINT 

NSTAR 

OD 

P 

P3 

POIS 
PINT 
PVAL 


IDIMX  +  1 

counter  used  in  storing  X3(I) 

integer  parameter  for  preparing  headings 

metal  thermal  conductivity  at  tempera- 
ture Z 

coefficients  used  in  determining  K(Z) 

K (X (I) ) *A(I) ,  useful  combination 

degree  of  polynomial  desired 

integer  parameter  for  preparing  headings 

LORDER  +  1 

counter  used  in  storing  significant 
stress  for  each  time  interval 

counter  used  in  storing  significant 
stress  of  adjoining  pipe  for  each  time 
interval 

input  from  data  card 

number  of  points  used  to  make  up  poly- 
nomial for  temperature  distribution 

1  (in  the  present  application  of  BETTY1) 

used  to  distinguish  two  joining  pipes 

number  of  time  intervals  between  lines 
of  printed  output 

option  parameter  in  data  deck 

outside  diameter  of  pipe 

a  useful  parameter  which  takes  care  of 
conversion  units 

quantity  useful  in  obtaining  end  moment 
and  radial  force 

Poisson's  ratio 

subroutine  used  to  integrate  polynomial 

subroutine  used  to  evaluate  /Trdr 
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Q 

R 

RAD (I) 

ROCK 

RT 
RY 

S3 

SHEAR 

SIG(I) 

SIGMAR 
SIGMAT 

SIGMAZ 
SIGTX 

SIGZX 

SIGR 

SIGT 

SIGZ 

SPHT(Z) 

SPHA,  SPHB, 
SPHC,  SPHD 

SSIG 

ST1 


quantity  useful  in  calculating  end 
moment  and  radial  shear  force 

radial  coordinate  of  pipe 

variable  used  to  store  radius  of  sig- 
nificant stress 

radius  of  maximum  significant  stress  for 
a  given  time 

absolute  value  of  ST1  -  ST2 

radial  distance  measured  from  mid-surface 
of  pipe 

quantity  useful  in  calculating  end 
moment  and  radial  force 

radial  force  as  a  function  of  axial 
coordinate 

variable  used  to  store  maximum  signifi- 
cant stress  for  each  time  interval 

radial  stress  remote  from  pipe  junction 

circumferential  stress  remote  from  pipe 
junction 

axial  stress  remote  from  junction 

circumferential  stress  due  to  bending 
moment  and  radial  force 

axial  stress  due  to  bending  moment  and 
radial  force 

total  radial  stress 

total  circumferential  stress 

total  axial  stress 

metal  specific  heat  at  temperature  Z 

coefficients  for  calculating  SPHT(Z) 


maximum  significant  stress  encountered 
principal  stress  1 
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ST2 
ST3 
STRESS (I) 

SURFCO(Z) 


SURA, 

SURB, 

SURC, 

SYRD 

T 

TAU 

TIMA 

TIME  (I) 

TMAX 

TOM 

TT1 

TT2 

TT3 

TVAL 
U 

UAY(I) 

UU 

UUU 

V 

VO 

w 


principal  stress  2 

principal  stress  3 

variable  used  to  store  significant 
stress  for  a  given  radius 

inner  surface  film  coefficient  at 
temperature  Z 

coefficients  for  calculating  SURFCO(Z) 

problem-time 

shearing  stress 

time  at  which  maximum  significant  stress 
in  encountered 

variable  used  to  store  time  of  maximum 
significant  stress  for  each  time  interval 

problem-end- time 

value  of  DT  for  first  pipe 

time  for  which  value  of  FELTA1  is 
recorded 

time  for  which  value  of  FELTA2  is 
recorded 

temperature  as  obtained  from  polynomial 
for  temperature  distribution 

subroutine  used  to  evaluate  TT3 

radial  displacement  remote  from  pipe 
junction 

alphabetic  data  for  headings 

value  of  right  hand  side  of  equation  (39) 

value  of  right  hand  side  of  equation  (40) 

deflection  at  pipe  junction 

radial  force  at  pipe  junction 

quantity  useful  in  calculating  end 
moment  and  radial  force 
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WT  wall  thickness 

X(I)  temperature  at  center  I 

XAVG  average  temperature  throughout  thick- 

ness of  pipe 

XB  value  of  z  where  bending  moment  is 

extremal 

X2(I)  data  point  denoting  radius  used  in 

fitting  a  polynomial 

X3(I)  coefficient  for  polynomial  for  tem- 

perature distribution 

XFO  outside  fluid  temperature 

XX (I)  newly  calculated  value  of  X(I) 

Y(I)  coefficients  of  polynomial  for  tem- 

perature after  integration 

Z  a  useful  quantity  for  calculating  XAVG 

ZED  used  in  FLUTEM 

ZR  absolute  value  of  ST3  -  ST1 

Recall  that  there  are  also  a  number  of  other  variables  with 
names  similar  to  those  above  but  with  a  terminal  1. 
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13  ABSTRACT 


A  method  is  devised  for  the  transient-state  thermal 
stress  analysis  of  two  pipes,  joined  butt  to  butt  and 
subjected  to  rapid  or  sudden  change  of  internal  fluid 
temperature.   Although  it  is  assumed  that  there  is  sym- 
metry about  the  common  pipe  axis,  the  properties  of  the 
materials  as  well  as  the  thickness  of  each  pipe  may 
be  different.   Then,  given  a  specified  time-temperature 
relationship  for  the  internal  fluid,  over  a  specified 
problem  time,  the  maximum  stress  encountered  may   be 
obtained.   A  digital  computer  program  is  appended  for 
the  solution  of  such  problems. 
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